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ABSTRACT 

Retrieval of orbital parameters of extrasolar planets poses considerable statistical chal- 
lenges. Due to sparse sampling, measurement errors, parameters degeneracy and mod- 
elling limitations there are no unique values of basic parameters such as period and 
eccentricity. Here we estimate the orbital parameters from radial velocity data in a 
Bayesian framework by utilising Markov Chain Monte Carlo (MCMC) simulations 
with the Metropolis-Hastings algorithm. We follow a methodology recently proposed 
by Gregory and Ford. Our implementation of MCMC is based on the object oriented 
approach outlined by Graves. We make our resulting code, ExoFit, publicly available 
with this paper. It can search for either one or two planets as illustrated on mock 
data. As an example we re-analysed the orbital solution of companions to HD 187085 
and HD 159868 from the published radial velocity data. We confirm the degeneracy 
reported for orbital parameters of the companion to HD 187085 and show that a low 
eccentricity orbit is more probable for this planet. For HD 159868 we obtained slightly 
different orbital solution and a relatively high 'noise' factor indicating the presence 
of an unaccounted signal in the radial velocity data. ExoFit is designed in such a 
way that it can be extended for a variety of probability models, including different 
Bayesian priors. 

Key words: stars:planetary systems - stars: individual: HD 187085, HD 159868 - 
techniques: radial velocities - methods: data analysis - methods: statistical - methods: 
numerical. 



1 INTRODUCTION 

The study of extra-solar planets is of great importance be- 
cause it helps us to understand the origin and evolution of 
the Solar System. The information gained from analysing 
new stellar systems serves as a sign post towards finding 
possible life forms outside our Solar System. The first 200 
or so detected extra-solar planets has been discovered us- 
ing the radial velocity method. Other detection methods in- 
clude astrometric method, transit photometry, gravitational 
micro-lensing and direct imaging. Astrometric method look 
for the periodic shifts in the position of star and transit 
method measures the attenuation of star light caused by the 
passage of the planet across the star. Gravitational micro- 
lensing utilises the amplification of light rays from the back- 
ground object by a an intervening massive object. In a num- 
ber of cases, existence of the planet has been confirmed by 
more than one method. The number of multi-planet systems 
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discovered has also increased during the past few years due 
the improved precision of radial velocity measurements. 



In this article, we are concerned with the extraction of 
orbital parameters from observed radial velocity data. A sur- 
vey of the literature in this field indicates that orbital param- 
eters and their uncertainties were traditionally obtained by 
a two step method. Searching for periodicity in the radial ve - 
locity data using a Lomb-Scargle (|Lomblll976l : IScarglelll982D 
periodogram to fix the orbital period and then estimat- 
ing o ther para meters by Levenberg-Marquardt (Leven berS 
1 19441 : iMarguardtl 1 19631 ) mini misation. Recent l y, Bayesian 
methods have been appfied (|Gregorvl l2005al : iFordI |2005| : 



iFord fc Gregonl l2007h to find out the best fit orbital pa- 
rameters and these studies suggest that Bayesian techniques 
have considerable advantages over the traditional methods; 
for e.g when the data does not cover a single orbital phase of 
the planet. In this work we analyse the radial velocity data 
from a Bayesian point of view and extract the orbital param- 
eters and corresponding uncertainties using Markov Chain 
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Monte Carlo (MCMC) simulations. We make our code, Ex- 
oFit, and a detailed documentation available on-lin^B- 

The outline of the paper is as follows. In Section [2] we 
summarise the modelling of radial velocities. In SectionOwe 
discuss the Bayesian approach to the problem. In Section |4] 
we present the MCMC implementation and the ExoFit soft- 
ware package. ExoFit is then applied to mock and observed 
data in Sections [5] and [6] respectively. The results are dis- 
cussed in Section [T] 



at rObservatoire de Haute Provence and La Silla Observa- 
tory (the Geneva extrasolar planet search). Radial velocity 
data for a star consists of time of observation ti, measured 
radial velocity Vi and uncertainty associated with each mea- 
surement Bi. These uncertainties are a characteristic of the 
instruments used for measurements. The precision of these 
instruments have im proved from the order of 10ms~^in 1 994 
to order of Ims"^ (jsutler et alll2006l : [Pepe et al.ll2004 ) at 
present This is extremely significant for finding low mass 
companions as well as planets with large a s. 



2 MODELLING OF RADIAL VELOCITY 



2.1 Doppler Spectrography 

Planets are many times fainter than their host stars be- 
cause they shine only by reflecting the star light. This makes 
their direct imaging extremely difficult. However, the grav- 
itational pull of the planet makes the star wobble and this 
produces measurable periodic shifts in the apparent speed 
of the parent star. The motion of the star around the cen- 
tre of mass causes the observed spectrum of the star to be 
Doppler shifted according its radial velocity, i.e. the velocity 
along the line of sight of the observer. This is measured over 
a course of time to obtain the radial velocity data along with 
the measurement uncertainties. 



2.2 Radial Velocity of Star 

A single planet model is assumed here to analyse the radial 
vel ocity data. The radial ve l ocity of a star ca n be written 
as l|Murrav fc Dermott|[200ol : lOhta et al.ll2005l ) 



V — sin(/i -I- zu) + e sin zu) 



(1) 



K 



where Vi is the ith radial velocity entry corresponding to 

time coordinate ti and, 

V = the systematic velocity of the system, 

mp = the mass of the planet, 

rUs = the mass of the star, 

n = ^ the mean motion and T is orbital period of planet, 

a = the length of the semi-major axis of the planet, 

i = the inclination of the orbital plane with the ecliptic, 

e = the eccentricity of the planet, 

fi = the true anomaly at time ti and 

zu = the longitude of periastron. 

The radial velocity depends on time via the true anomaly 
fi. The full formalism and comparison with a common ap- 
proximation is discussed in the User's Guide to ExoFit. 



2.3 Radial velocity data 



According to (http://exopIanet.eu) eighteen different radial 



velocity search programmes are looking for extrasolar plan- 
ets. Majority of the contributions come from Keck, Lick and 
Anglo- Australian observatories (the California & Carnegie 
and Anglo-Australian planet searches) and searches based 



3 BAYESIAN RETRIEVAL OF ORBITAL 
PARAMETERS 

3.1 Introduction 

The extraction of orbital parameters from the radial veloc- 
ity data poses considerable statistical challenges. Earlier in 
this article we mentioned the traditional two step method 
that is generally used to retrieve orbital pa rameters. Studies 
by [Cum ming. (2004 ) and I Gumming et al.l (|l9991 ) have iden- 
tified two cases where these methods become inefficient in 
accurately characterising the orbital elements: 

(i) When the orbital period is extremely short and the 
eccentricity is high. 

(ii) When the duration of observation does not span at 
least a single orbital phase. 

Incomplete radial velocity data gives rise to a multitude 
of orbital solutions which is referred to as parameter degen- 
eracy. Since the transit probability of the planet increases for 
short periods, the orbital parameters predicted by the pe- 
riodogram method can be verified, in some cases, with the 
help of transit photometry. Higher e ccentricities m a ke the 
radial velocity curve less sinusoidal. lO'Toole et all (|2007D 
makes use of a 2DKLS periodogram to incorporate the ef- 
fect of eccentricity of the orbits while searching for orbital 
per iods. Rec e ntly B aye sian t e chniq ues hav e been employed 
by iGregorvl (|2005al ). iFordI (|2005l ') and iFord fc Gregor^ 
l|2007l ) to retrieve the orbital parameters of extra-solar plan- 
ets. The results show that Bayesian methods tackle the dif- 
ficulties associated with the traditional methods efficiently 
and transparently. We emphasise that the relative merit 
of different methods depends on the quality of the data. 
Broadly speaking the choice of prior distribution may change 
the inference (posterior distribution) significantly when the 
quality of data is poor. 



3.2 The Bayesian method 

The starting poin t of any Bayesian analysis is Bayes' theo- 
rem fBavesl llTGSl '). Let y = (j/i, . . . , j/i, . . . , y„) be a vector 
of n observations whose probability distribution p(y|0, H) is 
conditional on k parameters = {6i, . . . ,0i, . . . ,9k), where 
H represents the background information or the hypothesis 
by which the probability statements are made. Suppose that 



^ http:/ /www. star. ucl.ac.uk/~lahav/exofit. html 



^ The measurement method and its uncertainties are discussed 
in the corresponding planet discovery papers. 
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the parameter 6 has the probabihty distribution p{9\H) 
Then, Bayes' theorem says 

9iy\e,H)p{e\H) 



p{0\y,H) = 



For a continuous 6 we can write 



p{y\H) = j p{y\e,H)p{e\H) dG , 



(2) 



(3) 



which is constant for given y and a probability distribution 
p{9\H). Then Equation [2] can be rewritten as 



p{0\y,H) = Cp{y\e,H)p{0\H). 



(4) 



In the above equation p{d\H) is called prior distribution of 
since it conveys our knowledge about 6 before the data 
has been observed. Correspondingly, p{d\y,H) is known as 
the posterior distribution of G given y. The factor C is a 
normalising constant which ensures that the posterior dis- 
tribution integrates to one. We call p{y\0, H) the likelihood 
function of since p(y\9, H) can be considered as a function 
of instead of y. Then, 



p{e\y,H) cc p(y\e,H)p{e\H). 



(5) 



Statistical inferences regarding 9 are derived from the pos- 
terior distribution of 9. The posterior distribution encapsu- 
lates all information about unknown quantities following 
the observation of the data y. The principal steps in the 
Bayesian method can be involve the Likelihood p{y\0,H), 
the prior p(9\ H), the posterior p(6 \ y,H) , and the resulting 
inference, e.g. lO'Hagan fc Forsteij l|2004 ). 



3.3 Likelihood function 

Let di represent the measured radial velocity data for the 
ith instant of time ti. Obs erved radial vel ocity data can be 
modelled by the equation (|Gregorvll2005ah 



di 



(6) 



where Ui is the true radial velocity of the star and is 
the uncertainty component arising from accountable but un- 
equal measurement errors which are assumed to be normally 
distributed. The term & explains any unknown measurement 
errors. There can be mu ltiple reasons for t he presence of this 
uncertainty component (|Butler et alll200^ V For example this 
could be the result of another planet in the system or caused 
by the intrinsic anomalies in the star spe ctrum due to the 
irregularities on the surface of the s tar l|Pepe et al. 1 12004 
iMavor et allbOOi iBouchv et al.ll2005h . Thus any noise com- 
ponent that cannot be modelled is described by the term 5. 
The probability distribution of 5 is chosen to be a Gaussian 
distribution with finite variance s^. Therefore the combina- 
tion of uncertainties ei + 5 has a Gaussian distribution with 
a variance equal to + s^. 

The true radial velocity Vi is modelled by the equa- 
tion ([2]). Thus the radial velocity Vi predicted by the math- 
ematical model at an instant ti is 

Vi = V — K sin(/i + w) + e sin -ro) . 

Six model parameters namely T, K, V, e, w,and X: as defined 
in the section (|2.2p are used to fit the above equation onto 
a given radial velocity data. 

Each error term in equation (|6]) is independent. Since 



Table 1. The assumed prior distribution of various parameters 
and their boundaries. It is similar to choice of priors given by 
iFord fc Greeorvl i2Q0l[) . except for the prior distribution of K. 

Para. Prior Mathematical Form Min Max 



T[days) 



Jeffreys 



Kims ) Mod. Jeffreys 

V ims~^) Uniform 

e Uniform 

■oj Uniform 

X Uniform 

s(ms~"'") Mod. Jeffreys 



-2000 2000 



1 

27r 
1 



In (in^ 



they are assumed to follow a Gaussian dis tribution, the lik e- 
lihood function is product of A'^ Gaussians (lGregorvl2005bl lah 
where A'^ is the number of observations. Thus 



P(y\0) = A exp 



E 



ii - Vi)"^ 



-2(0 



where 



A = (27r) 



-N/2 



2\-l/2 



(7) 



(8) 



and s becomes the seventh parameter in our probability 
model. 



3.4 Choice of Priors 

The choice of priors is extremely important in the Bayesian 
analysis as senseless choice of priors can produce to mis- 
leading results. Physical and geometric conditions govern 
the selection of prior distributions for most of the parame- 
ters. Since 9 = (T, K, V, e, ru, x, s) the prior distribution in 
our problem can be written as 



p{9\H) = p{T\H)p{K\H)p{V\H) 

p{e\H)p{zu\H)p{x\H)p(s\H) 



(9) 



on the assumption that they are independent. We will dis- 
cuss how the above conditions are met for our choice of prior 
for each parameter in the next few sections. 

We follow the choice of priors as given by 
iFord fc Gregonl (|2007D . as summarized in Table [T] 
Obviously, part of the Bayesian framework is the ability to 
change priors and to check the sensitivity of the results to 
them. Our ExoFit package allows this freedom. 



3.5 Posterior Distribution 

Posterior distribution is obtained by applying the Bayes' 
theorem given by the equation Useful and interesting 
features of the posterior distribution should be identified 
before making summary statements. For example, the poste- 
rior distribution may be unimodal but asymmetric or it can 
be multi-modal with many probability peaks. Any summary 
statistic(e.g. mean, median and mode) can be expressed 
in terms of posterior expectations of 9 l|Gilks et al.l Il996l : 
lBergeJ[l980l ). 
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4 MCMC IMPLEMENTATION 

4.1 The MCMC approach 

Difficulty in evaluating the multi-dimensional integrals is an 
inherent inability of any Bayesian formulation. Many tech- 
niques have been developed in the last 25 years to deal with 
this problem. Simulation methods dominate this area and 
several computational algorithms were developed to numer- 
ically integrate the posterior distribution in order to find 
out th e mar g inal d istributions of each parameter. Accord- 
ing to iBertj l|2004l ) the abundance of computational power 
has produced a paradigm shift with respect to statistics: 
Computationally intensive but conceptually simple methods 
are preferred. Markov Chain Monte Carlo (MCMC) method 
is one of the most commonly used methods for simulating 
complex probability distribution^. Our code is based on the 
concepts outlined bv iGraves. (,2007i ) . The emphasis is on the 
extensibility of the code to accommodate different probabil- 
ity models to a certain extent. 

Bayesian MCMC methods have gained popularity 
in various areas of astrophysics, for example in multi- 
paramet er estimation from cos mological data sets (e.g. Cos- 
moMC; iLewis fc Bridle! (|2002l )'). We emphasise that the two 
ingredients, the Bayesian approach and the MCMC tool, are 
distinct and not necessarily related. One may apply MCMC 
to a multi-parameter likelihood analysis (i.e. without priors- 
equivalent to a Bayesian method for uniform priors). On the 
other hand one may work out a full Bayesian method (with 
complicated priors) using a grid-based maximisation proce- 
dure, without the need for MCMC. Having said that, the 
combination of the Bayesian and MCMC methods is a very 
powerful way to tackle our problem. 

From a Bayesian point of view analysis of statistical 
problems requires an efficient tool for simulating posterior 
densities and MCMC methods are ideally suited for this 
purpose. In general, one must consider planet-planet inter- 
action while modelling radial velocity of a star. In typical 
systems the radial velocity of an n-planet model can be ap- 
proximated as a linear combination of n single planet radial 
velocities. With the present version of ExoFit it is possible 
to search for either one or two planets. 

4.2 The ExoFit software package 

ExoFit is a step towards achieving the goals mentioned 
above. It should be considered as a platform to develop 
MCMC based methods for estimating orbital parameters 
of a generalised multi-planet model. Object oriented design 
of ExoFit makes it extremely well suited for extending the 
analysis to multi-planet systems with prior constraints on 
several orbital parameters s uch as eccentr icity and length 
of semi-major axis. Following ICrave 3 (|2007D . our implemen- 
tation MCMC consists of Data, State, Bond and Update. 
They are referred to as objects in object oriented analysis. 
data handles the input data into the MCMC analysis. A 
state consists of a set of parameters whose posterior distri- 
bution is sought. The parameter values at a particular in- 

^ A gene ral form this meth o d given by Metropolis- 
Hastings llMetropoIis et al.l Il953l : iHastined Il970l) algorithm 
is explained in ExoFit User Guide. 



stant defines the state of Markov Chain in the analysis. The 
parameters in a particular state are connected to each other 
by a bond. It consists of prior densities and likelihood. For 
each state there corresponds a bond strength which is equal 
to prior x likelihood. In other words it is the posterior den- 
sity without the normalisation constant in Bayes' theorem. 
An update selects the parameters that should be updated 
at particular iteration. New values for the parameters are 
proposed according to the update defined and the new bond 
strength is then calculated for the proposed state. The new 
state is accepted or rejected according Metropolis-Hastings 
method. 

The central concept of this approach is that, the MCMC 
engine remains the same and need not be re-implemented 
whenever the probability model gets changed. We also take 
advantage of the commonalities among the different com- 
ponents of MCMC. Our implementation works for variety 
or prior distributions (Bonds) and Update methods. The 
only component that requires to be changed is the likeli- 
hood function. 

Convergence is an important aspect of any MCMC 
method and the choice of proposal distribution is absolutely 
crucial for achieving convergence of the chain. Choosing an 
effective proposal distribution is very difficult in many cases. 
E xoFit employs an ad aptive Metropolis algorithm described 
bv lHaario et al.l l| 19991 ) to fine tune the proposal distribution 
step sizes. 



5 APPLICATION TO MOCK DATA 

We used simulated data sets to test the accuracy of the or- 
bital parameters extracted by ExoFit. The output of ExoFit 
is analysed with the help of the R statistical environmenlQ 
which is freely available. A number of packages are avail- 
able in R for checking the convergence of Markov Chains 
and creating the histograms and density plots of posterior 
distribution of parameters. 

5.1 1-planet data with 1-planet fit 

In this The first set of radial velocity entries as shown in 
Table [2] (columns 1, 2 and 3) is created using a single planet 
model (equation (I2J) and adding a Gaussian noise. Table [3] 
shows the input parameters used to create the data set and 
the parameters obtained with ExoFit for a single planet 
search. Three posterior summary statistics are compared in 
Table O where column 2 gives the assumed input values, 
column 3 contains the mean values of MCMC samples for 
each parameter along with the sample standard deviation, 
column 4 sho ws the median val ues with 75% and 25% quan- 
tiles. We used iHvndmanI l|l996l) 's method to calculate High- 
est Density Regions (HDRs) of the posterior distribution. 
Column 5 contains posterior modes (maximum a posteri- 
ori) and associated 68.3% credible regions. The difference 
between actual radial velocity curve and the radial velocity 
curve for the parameters obtained from ExoFit is shown in 
the Fig. [1] We also plot the deviation of the radial veloci- 
ties from the actual one for each of these orbital solutions in 

* http://www.r-project.org/ 
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Table 2. Column 2 shows the simulated radial velocities obtained from a single planet model (Table [S]! using equation ((2J- A Gaussian 
noise with n = O.Oms"^ and a = 2.0ms~^ has been added to each radial velocity term. Columns 3 is obtained by adding the radial 
velocity from a second planet (see Table (Sjl to each term in column 2, thus creating a mock radial velocity data for a two planet system. 
The measurement uncertainties (Gaussian distribution with = 2.0ms~^ and a = 0.5ms~^) are shown in columns 4 and 5. 







J—J 1 1 \J 1 \ 1 1 LJ j 


IhV \l 1 to J 






1-planct 




2-planets 




1400.15200 


-17.90 


2.1 


-25.30 


2.1 


1413.87749 


-18.50 


2.8 


-18.20 


2.8 


1434.46573 


-30.00 


1.4 


-22.50 


1.4 


1455.05396 


-33.00 


2.0 


-26.20 


2.0 


1489.36769 


-41.10 


2.1 


-52.70 


2.1 


1496.23043 


-41.30 


2.4 


-50.80 


2.4 


1503.09318 


-45.50 


1.7 


-51.20 


1.7 


1564.85788 


-47.60 


1.7 


-45.30 


1.7 


1571.72063 


-44.20 


2.0 


-46.80 


2.0 


1585.44612 


-39.40 


1.6 


-50.70 


1.6 


1592.30886 


-33.60 


1.9 


-44.70 


1.9 


1612.89710 


-14.50 


1.9 


-14.70 


1.9 


1619.75984 


-6.17 


1.9 


-3.07 


1.9 


1654.07357 


38.20 


2.4 


45.30 


2.4 


1688.38729 


61.90 


1.6 


50.30 


1.6 


1695.25004 


64.90 


1.7 


54.90 


1.7 


1715.83827 


74.00 


2.9 


75.30 


2.9 


1722.70102 


73.70 


2.7 


78.00 


2.7 


1736.42651 


70.20 


1.8 


78.00 


1.8 


1743.28925 


72.70 


2.6 


81.10 


2.6 


1770.74024 


67.30 


2.3 


65.50 


2.3 


1777.60298 


65.10 


2.0 


57.80 


2.0 


1784.46573 


66.20 


2.7 


55.20 


2.7 


1798.19122 


61.90 


2.4 


53.40 


2.4 


1811.91671 


56.70 


1.9 


56.00 


1.9 


1818.77945 


56.20 


2.2 


58.80 


2.2 


1832.50494 


51.00 


1.7 


58.20 


1.7 


1901.13239 


32.80 


1.6 


25.90 


1.6 


1907.99514 


32.40 


2.0 


29.60 


2.0 


1921.72063 


29.80 


2.2 


33.70 


2.2 


1928.58337 


33.90 


3.7 


40.10 


3.7 


1935.44612 


29.10 


2.9 


36.80 


2.9 


1962.89710 


19.00 


2.2 


22.40 


2.2 


1969.75984 


17.40 


2.2 


16.30 


2.2 



Table 3. A comparison between input parameters to the mock data for 1 and 2 planets and the parameters extracted with ExoFit 
1-planet model. Columns 3, 4 and 5 show the posterior mean (and standard deviation), median (and 25% and 75% quantiles) and the 
maximum a posteriori, i.e. the posterior mode (and 68.3% highest density regions) obtained by the application of ExoFit 1-planet model 
to the single planet data from Table [2] Column 6 shows the orbital parameters of the second planet used to simulate the data for the two 
planet system (columns 4 and 5 in Table [2J. Summary statistics of the posterior distribution generated with ExoFit 1-planet model for 
this 2-planet data is shown in columns 7, 8 and 9. Note the difference between the noise factors extracted with ExoFit from the single 
planet data and the 2-planet data. Thus the noise factor s servers as an indicator for a possible embedded signal in the radial velocity 
data apart from the random Gaussian errors. 



Parameter 


Planot 1 


ExoFtt 
Mean 


ExoFtt 
Median 


ExoFit 
Mode 


Planet 2 


ExoFit 
Mean 


ExoFit 
Median 


ExoFit 
Mode 




12.0 


11.91 ± 0.44 


11 91+0-30 


11 91 + 0-44 




11.50 ± 1.60 


11 55+1-00 


11 69+l-4^ 


T(days) 


700.00 


704.85 ± 12.91 




702.79tii'6? 


100.0 


719.11 ± 53.06 


711.23+^1^= 


704 95+44'™ 
'^^■^^-43.50 


if (ms^l) 


60.00 


60.39 ± 0.56 


60 ^g+0.38 


60 "Sy+O-SS 


10.0 


62.52 ± 2.25 




62.48_2 2Q 


e 


0.38 


0.38 ± 0.03 


„ oo + O.Ol 
"■■^'-O-Ol 


„ „„ + 0.01 
"■•^°-0.02 


0.18 


0.40 ± 0.04 


40+" "^ 


40+'' "4 


zo 


3.10 


3.09 ± 0.03 


■5 nQ + 0-02 
•='■"^-0.02 


3 04 + " 


1.10 


2.96 ± 0.11 


2 96+°"^ 


2.96iS;ii 


X 


0.67 


0.67 ± 0.00 




67+°-'"' 
"■°'-o.oi 


0.17 


0.69 ± 0.01 




0.68tS;Sl 


s 




0.4276 ± 0.35 




09+0-43 
"■"^-0.09 




7.49 ± 1.14 


7 36+0-**0 


7 ,, + 0.93 
' --^^-0.93 
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Figure 1. Radial velocity curves created with mean, median and 
mode of the posterior distribution of parameters. The radial ve- 
locity curve for actual input parameters with the simulated data 
is also shown. 




1400 1500 1600 1700 1800 1900 2000 

Time (days) 



Figure 2. Deviation from the actual radial velocities of the sim- 
ulated data for 3 different summary statistics shown in Table [S] 
Note the differences are no more than 3ms~^. 

Fig. [2] We note the deviations between the estimators are 
smaUer than 3ms~^. The noise factor s « 0.4 is less than the 
measurement uncertainties, implying the assumed radial ve- 
locity model is "good" and there is no sign of an additional 
signal present in the data. 



5.2 2-planet data with 1-planet fit 

A second set of radial velocity entries as shown in column 
4 of Table [2] has been created by adding the radial velocity 
from a second planet with orbital parameters as shown in 
Table[3]to the radial velocity data of the single planet model 
in column 2 of Tabled The results from the application of 
ExoFit to this simulated data is shown in columns 7, 8 and 
9 of Table |3l Noise factor s ~ 7ms~^ in this case is clearly 
higher than the measurement uncertainties, indicating the 
presence of a signal unaccounted for in the radial velocity 
model. Note that the noise factor obtained in the previous 
case was only ~ 0.4ms~^. Fig. [3] shows the radial velocity 
curve for the median values of the posterior distribution of 
orbital parameters from ExoFit along with the actual radial 
velocity plot for the 2-planet model. 



5.3 A 2-planet fit to 2-planet mocli data 

ExoFit has an option to search for two planets in the ra- 
dial velocity data. We choose the probability model to be 
similar to that of the single planet model explained in Sec- 
tion |3]3l The observed radial velocity data is again modelled 
by the equation |6] The radial velocity Vi predicted by the 
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Figure 3. The radial velocity curve for two planet model with or- 
bital parameters of the planets from Table[3l along with the radial 
velocity plot created with median values of parameters obtained 
by ExoFit (median values) and the actual radial velocity. 



Table 4. The assumed prior distribution of orbital parameters 
and their boundaries for a 2-planet model. The boundaries for 
K2 can be made smaller in order to speed up the convergence of 
the Markov Chain. 



Para. 


Prior 


Mathematical Form 


Mill 


Max 


T^idays) 


Uniform 
Jeffreys 


1 


-2000 
0.2 


2000 
15000 


Xl(ms~^) 

ei 

XI 

T^idays) 


Mod. Jeffreys 
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Uniform 
Uniform 
Jeffreys 


(Kl-t-Jfio)~l 
... ( ^1 + ^1 max ^ 
> 

1 

1 

27V 
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1 

27r 
1 
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Uniform 
Uniform 

Mod. Jeffreys 


K2Q > 

1 

1 
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1 
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2000 
1 

27r 
1 

2000 



so 



mathematical model at an instant ti is given by 

Vi — Ki[sm{fii + zui) + ei sintui) 

+ K2{sin{fi2 + 1^2) + e2 smzu2)^ . (10) 

11 parameters {V,Ti, Ki,ei,wi,xi,T2, K2,e2,W2,X2'^ as 
defined in the Section [2.21 are used to fit the above equation 
onto the radial velocity data. The likelihood function is again 
given by equations [7] and [8] respectively, and s becomes the 
12**^ parameter in our probability model. The choice of prior 
distributions for each of these parameters is given in Table 

H 

We apply the ExoFit to the 2-planet mock radial veloc- 
ity data given in Table [2] (column 4). A comparison of actual 
orbital parameters with the ones extracted with ExoFit is 
given in Table 15.31 Fig 3] shows the actual radial velocity 
curve and the one created with the median values of the 
posterior distribution of orbital parameters obtained from 

^ Subscripts 1 and 2 indicate planets 1 and 2 respectively. 
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Table 5 . A comparison of actual values of orbital parameters used 
to create the mock radial velocity data for the 2-planet system 
(Columns 4 and 5 in Table [2J and the summary of the posterior 
distribution of orbital parameters extracted with ExoFit 2-planet 
model. Column 3, 4 and 5 shows posterior mean (and standard 
deviation), median (and 25% and 75% quantiles) and the maxi- 
mum a posteriori, i.e. posterior mode (and 68.3% highest density 
regions) respectively. 



Parameter 


Actual 


ExoFit 
Mean 


ExoFit 
Median 


ExoFit 
Mode 




12.0 


11.80 ± 00.52 






Ti(days) 


700.0 


709.06 ± 15.03 




708.03llt.36 


Ki (ms^ ^ ) 


60.0 


60.34 ± 0.62 


60.33l»fi 




ei 


0.38 


0.38 ± 0.01 




"•3S_o 01 


^1 


3.10 


3.10 ± 0.04 




, ,,+0.03 
•^■^^-0.04 


XI 


0.67 


0.67 ± 0.00 


-0.01 


67+°"° 


T2(days) 


100.0 


100.44 ± 0.53 


100.45l0-|5 


100.43±°j^ 




10.0 


10.19 ± 0.62 


in 1K + 0-.38 
10.18„^j 


10 14+°°= 




0.18 


0.19 ± 0.05 




20+° °'^ 




1.10 


1.27 ± 0.36 


, ,0 + 0-23 


1 29+°-^** 


X2 


0, 17 


0.15 ± 0.05 


is+°-03 

"■^''-0-03 


15 + °°'S 


s 




0.50 ± 0.41 


40+°-^^ 


15 + °°'S 
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Figure 5. Posterior density of 12 orbital parameters of the 2- 
planet model obtained from ExoFit. 



Figure 4. A comparison of the actual 2-planet radial velocity 
curve and the one obtained with ExoFit (median values). 



ExoFiE .Note that the noise factor s is now ^ 0.4ms 
compared to s ~ 7ms~^ in the single planet model. The 
density plots are shown in Fig [5] 

Constraining the longitude of periastron becomes in- 
creasingly difficu lt when the eccentricity of planetary system 
approaches zero. iFordI (|2006l ) suggests that, the efficiency of 
MCMC sampler in this situation can be increased by adopt- 
ing a new parametrisation with esin and ecos ro. For 
discussions about further improvements in parametrisation, 
choice of priors and proposal distributions for MCMC, see 



° We caution that, using median values for orbital parameters 
with skewed marginal posterior distributions can lead to models 
that very poor. 



iFordI (|2006l ). These considerations will be implemented in 
the future version of ExoFit. 



6 APPLICATION TO OBSERVED DATA 

In this section we apply ExoFit to published radial velocity 
data. We choose data sets in which the measurement noise 
is relatively high and the entries are poorly sampled. The 
following examples shows that ExoFit does a good job in 
estimating the posterior distribution orbital parameters. 



6.1 The companion to HD 187085 

Detection of a planet orbiting HD 187085 was announced in 
2006 with 40 epoch of radial veloci ty measurements b etween 
1998 November and 2005 October (| Jones et al.ll2006l ). It has 
reported that although an orbital solution of period 986 days 
and eccentricity 0.47 gave the best fit, a solution with low 
eccentricity also produced similar fit to the observed data. 
Application of ExoFit revealed the posterior density distri- 
bution orbital parameters as shown in Fig. |6] The MCMC 
values with iterations are shown in the user's manual. It is 
evident from Fig. |5] that the posterior density of eccentric- 
ity is a heavily skewed . The joint probability distributing 
of orbital period and eccentricity is shown in Fig. [7] with 
two dimensional HDRs. From these plots it is reasonable 
to conclude that the data clearly favours a low eccentricity 
orbit. 
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Figure 6. Posterior distribution of orbital parameters of the com- 
panion to HD 187085. These plots were created with R from the 
output of ExoFit. 



Table [6] shows the summary statistics of the posterior 
distribution of orbital parameters of HD 187085 along with 
the published results. The eccentricity is e = 0.33, 0.30, 0.11 
from the mean, median and mode respectively, compared 
with the published estimate e — 0.47. Keplerian orbital 
solutions for HD 187085 from Table |6] shown in Fig.[8l 

How sensitive are the results to the assumed priors? In 
Table [7] we show results based on 10 trial runs of for Gre- 
gory's priors (Table [l| and for Top Hat priors, i.e. uniform 
between the same lower and upper values as in Gregory's pri- 
ors. This indicates the data and assumed model are 'good'. 
The role of priors would be more dramatic in the case of 
poor data. 



6.2 The companion to HD 159868 

Orbital parameters of a planet orbiting HD 159868 from 
28 radial velocity m easurements was reported on 2007 by 
lO'Toole et al.l |20o3). They employed a two-dimensional pe- 
riodogram to search for the optimum orbital period and the 
eccentricity for the orbit. The posterior distribution of pa- 
rameters are shown in Fig. [9] and since they are nearly sym- 
metric, we choose median of the distribution as the point 



^ For a posterior distribution that is approximately Gaussian 
(e.g.. parameter T), estimates like mean, median and mode will 
yield nearly the same values. However, when the posterior density 
is not Gaussian and exhibit skewness (e.g.. parameter e) mean, 
median and mode will differ significantly. Since this is the case 
for parameter e, posterior median and and mode will be a better 
estimate than posterior mean. If the posterior distribution has 
more than one peak, posterior modes provide ideal summary of 
the distribution. 



Figure 7. The joint posterior distribution of orbital period T 
and eccentricity e for the companion to HD 187085. The con- 
tours represent probability coverage=0.01,0.05,0.1,0.2,0.3,0.4,0.5 
and 0.7. The circle in the centre represent the mode of the joint 
probability distribution. The dots outside the contour show the 
samples that fall outside 99% coverage probability. 



Median:T=1065.74 days, e=0.30 - 
Wode:T=1061.54days,e=0.11 - 
Radial yeldcity data ^ 




Figure 8. Possible orbital solutions for the companion to 
HD187085. 



statistic. Table |8]shows the median of samples from ExoFit 
and the Keplerian orbit obtained is plotted in Fig. 1 101 along 
with the observed data and uncertainties. Table[7]shows the 
sensitivity to changing the priors to Top Hat. We note that 
the 'noise' factor s is 10 times larger than the measurement 
errors, probably indicating the presence of an unaccounted 
signal in the radial velocity data. 



7 DISCUSSION 

We present a new software package, ExoFit, for estimat- 
ing the orbital parameters from radial velocity data irr a 
Bayesian framework by utilising Markov Chain Monte Carlo 
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Table 6. A comparison of estimators for orbital parameters of the companion to HD 187085. The columns provide the posterior mean 
(and standard deviation), median (and 25% and 75% quantiles) and the maximum a posteriori, i.e. posterior mode (and 68.3% credible 
regions). The last column is from the published results of Jones et al. (errors were not provided in their paper). Here rup sin i and a are 
calculated by assuming the mass of the star = 1.16 Mq. 



Parameter 


ExoFit 


ExoFit 


ExoFit 


Jones et al 




Mean 


Median 


Mode 




T{days) 


1065.84 ± 45.75 


1065.74l|°:^« 


IO6I.54+5I353 


986 


X(ms-l) 


16.55 ± 4.49 




15.Q8till 


17 


y(ms-i) 


-0.94 ±1.56 


-0.9311:°° 






e 


0.33 ± 0.22 


o.3ol°;S 


11+0-29 


0.47 


■UJ 


118.57 ±23.19 


112.39llfj8 




94 


X 


0.11 ±0.06 




10+007 




s 


5.47 ± 1.08 


5 40+0 ''2 






Tp (JD -2450000) 


993 ± 60.49 


992ti 




912 


rUp sin i(M J „p) 


0.81 ±0.10 


o.sitro? 




0.75 


a{AI ') 


2.11 = ().()() 




2.11';!:;!;; 


2.05 



Table 7. Mean and standard deviation of posterior modes of parameters obtained from 10 trial runs of ExoFit. 





HD 187085 


HD 187085 


HD159868 


HD159868 


Parameter 


Gregory's 


Top Hat 


Gregory's 


Top Hat 


T (days) 


1068.49 ± 3.94 


1068.57 ± 2.36 


998.04 ± 4.88 


1007.38 ± 5.23 


K (ms-i) 


15.35 ± 0.30 


16.34 ± 0.66 


43.16 ±4.93 


41.44 ± 2.00 


V (ms-i) 


-1.03 ±0.10 


-0.99 ±0.12 


1.13 ± 0.21 


1.05 ±0.28 


e 


0.14 ±0.03 


0.16 ±0.51 


0.70 ± 0.04 


0.68 ±0.02 


TO 


96.30 ±2.86 


96.87 ±2.29 


94.05 ± 4.28 


98.50 ±3.28 


X 


0.11 ± 0.004 


0.11 ±0.05 


0.68 ± 0.002 


0.68 ±0.001 


s {ms~^) 


5.28 ±0.07 


5.03 ±0.06 


9.52 ± 0.10 


9.77 ±0.15 


Tp (JD -2450000) 


992.96 ± 0.42 


993.74 ± 0.05 


706.15 ± 2.77 


689.35 ± 1.76 


rrip sin i {Mjup) 


0.80 ± 0.006 


0.81 ± 0.008 


1.61 ±0.03 


1.60 ±0.01 


a {AU) 


2.14 ±0.005 


2.14 ±0.003 


2.01 ± 0.006 


2.02 ± 0.006 



(MCMC) simulations with the Mctropohs-Hastings algo- 
rithm. ExoFit can search for one or two planets in a given 
radial velocity data and can easily be extended to search for 
more planets. Wo applied ExoFit to simulated data sets to 
check the accuracy of the parameters extracted. As an il- 
lustration we re-analysed the orbital solution of companions 
to HD 187085 and HD 159868 from the published radial ve- 
locity data. We confirm the degeneracy reported for orbital 
parameters of the companion to HD 187085 and show that 
a low eccentricity orbit is more probable for this planet. For 
HD 159868 we obtained slightly different orbital solution 
and a relatively high 'noise' factor indicating the presence 
of an unaccounted signal in the radial velocity data. We 
have also studied the sensitivity of the results to changes in 
the Bayesian priors. We plan to extend ExoFit to search for 



more than two planets and to analyse transit photometry 
data. 
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Table 8. A comparison of the orbital parameters of the companion to HD 159868 from ExoFit and the published results of[0 'Toole et al.l 
Here rrip sin i and a are calculated by assuming the mass of the star rris = 1.09 Mq. 



Parameter 


ExoFit 


ExoFit 


ExoFit 


O'Toole et al 




Mean 


Median 


Mode 




T (days) 


1012.28 ±34.53 


1008.25124.6^ 


997.05l^5;« 


986 ±9 


K (ms-i) 


46.87 ±25.63 


4i-30l|;J^ 


40.34l«;f 


43.3 ±2 


V (ms-i) 


0.96 ±2.80 




0.8912-71 
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0.67 ±0.11 
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Figure 9. Posterior distribution of orbital parameters of the com- 
panion to HD 159868. 



Figure 10. Best fit Keplerian orbit for the companion to HD 
159868. 
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